/************************************************************************************
**		tree.cpp -- FAST MULTIPOLE METHOD tree routines..
wellseperate changes...

findroot() fixed.

added createRepo()	initialize the Repository with 3.BlockSize
added mynew()	use a Node array..
***********************************************************************************/

#include "tree.h"
#include "multipole.h"
#include "mpole_calc.h"
#include "mpole_misc.h"


/************************************************************************************
**		clearNode, clear the content of the Node 
************************************************************************************/
int tree::clearNode ( Node* oldNode )
{
	int k;

/************* time ****************************************************************************/
#ifdef xdebug    		
	if ( oldNode == NULL ) { cout<<"clearNode icin null root ";return -1;}
#endif    
/*****************************************************************************************/
	
	for ( k=0; k < NSUB ; k ++ ) 
		oldNode->child[ k ] = NULL;
	
	oldNode->geocenter = 0;

	oldNode->index = 0;

	Clear_local( oldNode );				//oldNode->l;  local exansion
	
	oldNode->level = 0;

//	oldNode->lvalid = false;			//local valid flag 0 

	Clear_mpole( oldNode );				//oldNode->m multipole expansion
	
//	oldNode->mvalid = false;			// multipole valid flag
	
	oldNode->nneighbour = 0;			//number of near neighbour with self 2D=9,3D=27

	//for (i=0;i< NINTERACT; i++)				//NINTERACT = # of interaction list 
	//	oldNode->nearneighbour[ i ] = NULL;

	oldNode->nearneighbour.clear();
	
	oldNode->Nparticle = 0;

	//oldNode->plist();

	oldNode->rsize = 0;

	oldNode->Type = EMPTY;	

	oldNode->parent = NULL;
	
	return 0;
}

/************************************************************************************
**		clearRepo, allocate and initialize the Node's array.. 
**		it's a two dimensional array..
**		Nodelist[index1][index2]
************************************************************************************/
void tree::createRepo ( long N )
{
	Nparticle = N;	//number of particles.
	currentNode = 0;	//need for reuse of Nodelist every time step.
}

/************************************************************************************
**		mynew() return a free Node from the array.
**		
************************************************************************************/
Node* tree::mynew( void )
{
	if ( currentNode < (long)Nodelist.size() )
	{
		clearNode( &Nodelist[currentNode] );
		return ( &(Nodelist[currentNode++]) );
	}
	
	Node *returnthisNode = new Node;

	Calloc( &(returnthisNode->m),Mtermp);
	Calloc( &(returnthisNode->l),Mtermp);

	clearNode( returnthisNode );
	
	Nodelist.push_back(*returnthisNode);

	currentNode++;

	return returnthisNode;
}


/************************************************************************************
**		clearNodelist() reset the Node's 
**		
************************************************************************************/
void tree::clearNodelist()
{
	currentNode = 0;
}




/************************************************************************************
**		emptytree release the memory for all NODES 
**		free the tree memory 
************************************************************************************/
void tree::emptytree( char *msg )
{
sprintf(msg , "DirectForce;%d;M2M;%ld;M2L;%ld;L2L=%ld;timeM2L;%0.3f;timeL2L;%0.3f;timeNear;%0.3f;timeDirectForce;%0.3f;timeForceC;%0.3f;Addmultipole;%ld;"\
,countdirectforce, countM2M, countM2L, countL2L , timeM2L/(float)CLOCKS_PER_SEC,\
timeL2L/(float)CLOCKS_PER_SEC, timeNear/(float)CLOCKS_PER_SEC,timeDirectForce/(float)CLOCKS_PER_SEC,\
timeForceC/(float)CLOCKS_PER_SEC, countAddmultipole );

	countaddleaf = 0;
	countLeaf2Node = 0;
	countdirectforce   = 0;
	countAddmultipole = 0;
	countM2M = 0;
	countM2L = 0;
	countL2L = 0;

	clearNodelist();
	
	clearNode( FMMroot );

	FMMroot->index = 0;
	
	FMMroot->geocenter = ( 0 );

	FMMroot->rsize = 0;
}

//////////////////////////////////////////////////////////////////////
// TREE Construction
//////////////////////////////////////////////////////////////////////
tree::tree( particlePtr partList , double stepsize=0.01 , long N=1 ) 	
{	
	treepartList = partList;		//particle list array 
	dt=stepsize;					//stepsize of simulation needed for force calculation  ??? 

	nlevel = 0;						//number of level of tree

	Nparticle = N;
	
	FMMroot = new Node;				//root of the FMM tree = FMMroot
	Calloc( &FMMroot->m , Mtermp );	//FMMroot not need multipole and local
	Calloc( &FMMroot->l , Mtermp );
	clearNode ( FMMroot );

	//default values
	FMMroot->geocenter = 0;
	FMMroot->rsize = 0;
	//geocenter and rsize will be replaced by findroot function

	//tailNode=headNode=NULL;
	
	createRepo( Nparticle );

	countaddleaf = 0; 
	countLeaf2Node = 0;
	countdirectforce   = 0;
	countAddmultipole = 0;
	countM2M = 0;
	countM2L = 0;
	countL2L = 0;


	timeTemp = 0;
	timeM2L = 0;
	timeL2L = 0;
	timeDirectForce = 0;
	timeForceC = 0;
	timeNear = 0;

	for ( int i=0;i<Nparticle;i++)         //give all particle an ID
   {
      treepartList[i].Id = i;
   }
}


/************************************************************************************
**		load particles into tree , start from FMMroot node use the particle array partList
**		and the number of particles N .
************************************************************************************/
void tree::loadparticles( )
{
	long i;

	findroot( Nparticle , treepartList , FMMroot->rsize , FMMroot->geocenter );

	int j;
	bool isinserted = false;
	Node *root;

	//insert particle into the root cell one by one..

	for ( i=0; i<Nparticle; i++)
	{
		//start from root of the tree and traverse the tree.
		root = FMMroot;

		while ( isinserted == false )
		{
			j = addparticle( i , &treepartList[i] , root );
			if ( j == -1 ) isinserted = true;
			else if ( j>=0 && j< NSUB ) root = root->child[j];		//go one level down, child node
			else cout<<"return value of addparticle must be between -1 and 8..";

		}//while

		isinserted = false;

	}

	//displaytree( FMMroot );

	//cout<<endl<<"end of load particles.."<<endl;
}

/************************************************************************************
**		add particles into tree recursively , start from root node 
**		
************************************************************************************/
int tree::addparticle ( long index , particle *part , Node* root )		
{			
	//which child of the root, particle will place  ??
	int wcld = whichchild( part->getpos() , root );

	Node *childcell = root->child[ wcld ];

	if ( childcell == NULL || childcell->Type == EMPTY ) 
	{
		addleaf( index  , part , root , wcld );
		root->Nparticle++;

		//leaf's parent is root..
		root->child[wcld]->parent = root;

		return -1;//particle inserted.. added.
	}
	else if ( childcell->Type == LEAF )
	{
	
		//look for # of particle in LEAF
		if ( childcell->Nparticle < particleLEAF )
		{
			childcell->Nparticle++;				//number of particle in childcell
			childcell->plist.push_back(part);	//add a pointer to this particle

			part->parent = childcell;
			root->Nparticle++;
			return -1;		//particle inserted return ok.
		}
		else //# of particle in LEAF is bigger then NPART
		{
			leaf2node( root , wcld );			//convert root's wcld child to a node.
			root->Nparticle++;
			return wcld;
		}
	}

	else if ( childcell->Type == NODE )
	{
		root->Nparticle++;
		return wcld;
	}
	cout<<"there is a problem in addparticle";
	return -2;
}


/************************************************************************************
**		which child of the root NODE will place the particle with position pos 
**		and return the number of child node
************************************************************************************/
int tree::whichchild( vektor pos , Node* root )
{
	if ( root == NULL ) 
	{
		cout<<"in whichchild function root is a NULL value..";
		return -2;
	}

	int childnumber=0;

	for ( int k=0; k<NDIM; k++ )
	{
		if (pos[k] >= root->geocenter[ k ] )
		{
			childnumber += 1<<k;
/************* extra debugging ****************************************************************************/
//#ifdef xdebug 
//			if ( pos[k] > (root->geocenter[k] +	root->rsize / 2)   )
//			{
//				cout<<"sorun var .. "<<pos[k] << " root " << root->geocenter[k]<<" rsize="<<root->rsize<<" toplam = " << (root->geocenter[k] +	root->rsize / 2) <<endl;;
//				return -2;
//			}
//
//		}
//		else if ( pos[k] < ( root->geocenter[k] -	root->rsize / 2)   )
//		{
//			cout<<"pos root icinde degil.." << pos[k] << " root " << root->geocenter[k]<<endl;
//			return -2;
//#endif
/************* extra debugging ****************************************************************************/

		}
	}//for

	if ( childnumber>=0 && childnumber<8 ) 
		return childnumber;
	else 
	{
		cout<<"childnumber'da sorun var = "<<childnumber;
		return -2;
	}
	
	
}

/************************************************************************************
**		add leaf to the NODE root 
**		must to know which child of the root..
************************************************************************************/
void tree::addleaf( long index , particle* part , Node* root , int wcld)
{

/************* extra debugging ****************************************************************************/
#ifdef xdebug 
	// insert a leaf to root's wcld child the particle.	
	if ( root->Type == LEAF )
	{
		cout<<"\naddleaf LEAF node 'a ekleme yapiyor.";
		display( root );
	}
#endif
/************* extra debugging ****************************************************************************/

	root->Type = NODE;					//root was EMPTY and now it is a NODE

	Node* newLeaf = mynew();
	
	newLeaf->level = root->level+1;

	newLeaf->Nparticle++;
	newLeaf->plist.push_back(part);
	
	part->parent = newLeaf;

	newLeaf->rsize = root->rsize / 2.0;
/************* extra debugging ****************************************************************************/
#ifdef xdebug 
	if ( newLeaf->rsize == 0 ) 
		cout<<"add leaf icinde rsize hesaplamasinda 0 var.. L="<<newLeaf->level<<" index="<<newLeaf->index;
#endif
/************* extra debugging ****************************************************************************/
	
	newLeaf->Type = LEAF;

	/* calculate the center of the leaf cube..	*/
	for ( int k =0; k < NDIM ; k++ ){
         if ( ( wcld >>k) % 2 ) 
			 newLeaf->geocenter[k] = root->geocenter[k] + newLeaf->rsize/2.0;
         else newLeaf->geocenter[k] = root->geocenter[k] - newLeaf->rsize / 2.0;
	}
	
	/* newLeaf ok */ 
	
	root->child[wcld] = newLeaf;

	//number the lowest level of the tree..
	if ( nlevel < newLeaf->level ) nlevel = newLeaf->level; 

}

/************************************************************************************
**		convert leaf to a NODE  root->child[ wlcd ] is a leaf convert to node.
**		convert root nodes wcld'th number of children leaf to a node
**		add the old leaf to the new node as a leaf..
************************************************************************************/
void tree::leaf2node( Node* root , int wcld )		
{
	Node* oldleaf = root->child[ wcld ];
	
	countLeaf2Node++;

	int Nparticle2 = oldleaf->Nparticle;
	
	oldleaf->Nparticle = 0;

	vector<particle*> tempplist = oldleaf->plist;

	oldleaf->Type = NODE;

	vector<particle*>::const_iterator temppart; 

	for (temppart = oldleaf->plist.begin(); temppart != oldleaf->plist.end(); ++temppart) 
	{
		wcld = whichchild( (*temppart)->getpos() , oldleaf );
		addparticle( (*temppart)->getId() , *temppart , oldleaf );

	}

	oldleaf->plist.clear();

}

/*************************************************************************
**      Upward PASS of the tree 										** 
**		Calculate multipole expansion of particles and					**
**		shift multipole expansion to parent node 
**************************************************************************/
void tree::upwardpass( void )
{
	//if first pass, allocate space for variables.
	static bool firstpass = true;
	
	if ( firstpass == true )
		/* setup multipole expansion , allocation space for series */
		MultipoleSetup( Mtermp );

	/*calculate multipole expansion */ 

	//calc multipole expansion only for particles in LEAVES
	calMultipole();

	//calc multipole expansion by shifting exp.
	calMultipole2( FMMroot );

	firstpass = false;

	//cout<<endl<<"end of upward pass ..."<<endl;
	//exit(3);
	//c+out<<endl<<"tree "<<endl;
	//displaytree( FMMroot );

} 


/*************************************************************************
**		Bottom-up Traversal												**
**		Calculate multipole expansion of particles and					**
**		shift multipole expansion to parent NODE						** 
**	if root is a LEAF calculate ME of LEAF parent			 			**
**************************************************************************/
void tree::calMultipole() 
{
	long i;
	vektor vect=0;
	Node *parentLeaf;

	for ( i=0; i<Nparticle;i++)
	{
		//add i. particle multipole exp. to its parent LEAF
		parentLeaf = treepartList[i].parent;

		//if this child is a LEAF calculate multipole expansion with partcile in this LEAF
		if ( parentLeaf->Type == LEAF )
		{
			
			/* vec is the distance between the particle and center of root Cell*/
			vect = treepartList[i].getpos() - parentLeaf->geocenter;
			AddMultipoleC( parentLeaf->m , Mtermp , treepartList[i].getdata() , vect );
			countAddmultipole++;

			//Display_mpole( parentLeaf );
			//parentLeaf->mvalid = true;
		}//if
/************* extra debugging ****************************************************************************/
#ifdef xdebug 
		else 
		{
			cout<<"calcMultipole icinde addmultipole cisimleri eklerken parent LEAF degilmis.";
			display( parentLeaf );
		}
#endif
/************* extra debugging ****************************************************************************/


	}//for i.

}

/*************************************************************************
**		Calculate multipole expansion of particles and					**
**		shift multipole expansion to parent NODE						** 
**	if root is a LEAF calculate ME of LEAF parent			 			**
**************************************************************************/
void tree::calMultipole2( Node* parent ) 
{
	int i;
	Node* childNode;
	
	// for all children of parent NODE 
	for( i=0; i<NSUB; i++)				
	{
		// childNode is the child of parent
		childNode = parent->child[i];

		// always control for NULL
		if ( childNode != NULL  && childNode->Type!= EMPTY )
		{

			if ( childNode->Type == LEAF )
			{
				Calc_M2M( childNode , parent );
				countM2M++;
	
				//cout<<"calc_M2M if LEAF sonrasinda..."<<childNode->parent->index<<endl;
				//Display_mpole( parent );
				//MathdumpY_C( parent->m , Mtermp, "matdump.txt", 1 );
				
			}//if
			

			/*if its a NODE, than recursively calculate all childrens Multipole expansion
			and shift all children multipole expansion's to parent cell */
			else if ( childNode->Type == NODE )
			{
				calMultipole2( childNode );		//recursive calc.
				
				Calc_M2M( childNode , parent );
				countM2M++;
			}//else if
		}//if NULL
	}//for

}



/*************************************************************************
**      Downward PASS of the tree 										** 
**		Calculate local expansion for cells and							**
**		shift local expansion to parent node							**
**************************************************************************/
void tree::downwardpass( )
{


	int i;
	//downwardpass start from root cell and root has no local expansion 
	//root->children has no local expansion but have near neighbours 
	//so calc root->children near neighbour manually.
	buildNearRoot( );

    //root->children->children has local expansion 
	//calculate local expansion for all nodes and shift to children
	clock_t start = clock();

	for ( i=0; i<NSUB; i++ )
		if ( FMMroot->child[i] != NULL  && FMMroot->child[i]->Type!= EMPTY  ) 
			calcLocal( FMMroot->child[i] );
/************* time ****************************************************************************/
#ifdef showtime    		
	cout<<"\ntimeM2L="<<timeM2L / (float)CLOCKS_PER_SEC<< 
			" -- timeL2L = "  << timeL2L / (float) CLOCKS_PER_SEC << 
			" -- timeNear Neighbour = " << timeNear / (float)CLOCKS_PER_SEC<< endl;
	cout<<"\ncalculation of local expansion takes "<< ( clock()-start) / (float)CLOCKS_PER_SEC<<" seconds"<< endl;
#endif    
/*****************************************************************************************/
	

}


/************************************************************************************
**		Build near neighbour list so the interaction list may be calculated			**
**	this will calculate only for FMMtree's root and its children.
** FMMtree's children has only near neighbours , no interaction list.
************************************************************************************/
void tree::buildNearRoot( )
{
	int i=0,j=0;
/************* extra debugging ****************************************************************************/
#ifdef xdebug 
	if ( FMMroot == NULL ) { cout<<"big trouble..FMMtree's root NULL"; }
#endif
/************* extra debugging ****************************************************************************/

	Node *childNode,*otherChild;

	//for all children of root NODE 
	for( i=0; i<NSUB; i++)				
	{
		//childNode is the child of root
		childNode = FMMroot->child[i];

		// always control for NULL
		if ( childNode != NULL )
		{
			//calculate near neighbour of this child..
			for( j=0; j<NSUB; j++)
			{// always control for NULL
				otherChild = FMMroot->child[j];
				if ( otherChild != NULL )
				{
					//self is near neighbours so add itself.
					childNode->nneighbour++;		//number of neighbour
					childNode->nearneighbour.push_back(otherChild);
					
				} // if NULL
				else if ( childNode->Type== EMPTY ) cout<<"root'un cocuklarindan EMPTY olan var.";
			}//for j 

		} // if NULL
		//else if ( childNode->Type== EMPTY ) cout<<"root'un cocuklarindan EMPTY olan var.";
	}//for i 
}


/************************************************************************************
**		calc Local expansion recursively start from root 
************************************************************************************/
void tree::calcLocal( Node* parent )
{
	int i,j,k;
	 /*if its a NODE, than recursively calculate all childrens Multipole expansion
	 and shift all children multipole expansion's to parent cell */
	Node *childNode;

/************* extra debugging ****************************************************************************/
#ifdef xdebug 
	if ( parent == NULL ) cout<<"calcLocal icinde root NULL";
	//else if ( parent->Type == LEAF ) { cout<<"calcLocal has a LEAF parent";exit(12);return; }
#endif
/************* extra debugging ****************************************************************************/
	if ( parent->Type == NODE )
	{
		for ( i=0; i<NSUB; i++ )
		{
			//for all children of the root cell 
			 childNode  = parent->child[i];

			  if ( childNode != NULL  && childNode->Type!= EMPTY  )
			 {		
				 timeTemp = clock();

				 // look for parent cell's near neighbours for interaction list.
				 for ( j=0;j<parent->nneighbour; j++ )
				 {
					 interact( childNode , parent->nearneighbour[j] );
				 }//for j

				 timeNear += clock() - timeTemp;

				 timeTemp = clock();

				 //shift local expansion to all children
				 for ( k=0;k<NSUB;k++)
					 if ( childNode->child[k] != NULL && childNode->child[k]->Type!=EMPTY) 
					 {
						Calc_L2L( childNode , childNode->child[k] );
						countL2L++;
					 }

				timeL2L += clock() - timeTemp;
			 }//if NULL
		 }//for i

		 // calculate local expansion recursively.
		 for ( i=0; i < NSUB ; i++ )
		 {
			 childNode  = parent->child[i];

			 if ( childNode != NULL  && childNode->Type!= EMPTY && childNode->Type!=LEAF)
			 {
				 calcLocal( childNode );
			 }
		 }//for i

	 }//if NODE
}

/************************************************************************************
**	wellseperate ??
*************************************************************************************/
bool tree::wellseperate( Node thiscell , Node othercell )
{
	double r = dist ( thiscell.geocenter , othercell.geocenter );
    r = sqrt( r );
	double length = thiscell.rsize + othercell.rsize;

	//length = length * SQRT3;
	length = length * SQRT3 ;

	if ( r > length )
	   return true;
   else 
	   return false;
}


/************************************************************************************
**		interact this cell with the thiscell's parent's near neighbours
**	neparent = is the near neighbour of thiscell 's parent so look for all children of
**  neparent and interact with or add near neighbours
************************************************************************************/
void tree::interact( Node* thiscell , Node* parentNear )
{
	//interact this cell with thiscell's parent cell's near neighbours = parentNear 
	//thiscell and parentNear is not at same level of tree
/************* extra debugging ****************************************************************************/
#ifdef xdebug    		 
	if ( thiscell== NULL ) cout<<endl<<"in interact thiscell NULL";
	if ( parentNear== NULL ) cout<<endl<<"in interact neparent NULL";
#endif
/************* extra debugging ****************************************************************************/

	int i=0;
	Node *childNode;

	if ( wellseperate( *thiscell , *parentNear ) == true ) 
	{
		timeTemp = clock();

		Calc_M2L( parentNear , thiscell );
		countM2L++;
		//cout<<endl<<"interact thiscell with parentNear";
		//interaction list not neighbour
		timeM2L += clock() - timeTemp;
	}
	else		//interact with all children of neparent 


	{

		if ( parentNear->Type == LEAF ) 
		{
//////////////////////////////////////////////////////////////////////////////////////////////
//bu if sartina girmiyor..
			
			if ( wellseperate( *thiscell , *parentNear ) == true )
			{
				timeTemp = clock();
				//convert childNode multipole expansion to thiscell local expansion.
				Calc_M2L( parentNear , thiscell );
				timeM2L += clock() - timeTemp;
			}//if wellseperate
			else //not wellseperate
			{
			
				timeTemp = clock();

				thiscell->nneighbour++;
				thiscell->nearneighbour.push_back(parentNear);
				
				timeNear += clock() - timeTemp;
//////////////////////////////////////////////////////////////////////////////////////////////
			}//else

		}//if LEAF
		
		else if ( parentNear->Type == NODE )
		{
			/* for all children of neparent NODE */ 
			for( i=0; i<NSUB; i++)            
			{/* always control for NULL */
				childNode = parentNear->child[i];
				if ( childNode != NULL  && childNode->Type!= EMPTY )
				{/* newNode is the child of root */ 
					//for interact look for the distance between
					if ( wellseperate( *thiscell , *childNode ) == true )
					{
						timeTemp = clock();
						//cout<<endl<<"convert Multipole expansion of newNode to thiscell local exp.";
						//convert childNode multipole expansion to thiscell local expansion.
						Calc_M2L( childNode , thiscell );
						countM2L++;
						//cout<<"interact icinde =";
						//Display_local( newNode );
						timeM2L += clock() - timeTemp;
					}//if wellseperate
					else // not wellseperate
					{
						timeTemp = clock();

						thiscell->nneighbour++;
						thiscell->nearneighbour.push_back(childNode);

						timeNear += clock() - timeTemp;
					}//else
				}//if !NULL
			}//for i
		}//else if parent NODE

	}//else	interact with all children of neparent 
}

/************************************************************************************
**		for calculating forces, first we need to calculate local expansions,
**		than with interaction list top-down traversal tree and calculate children local
**		expansion.
**		
************************************************************************************/
void tree::calculateforces( )
{
	int k=0;
	long i;
	vektor vect=0;

	//control the ID of all particles.ID is the index of the particle array.
	for (i=0;i<Nparticle;i++)
	{
		treepartList[i].force = 0;
/************* extra debugging ****************************************************************************/
#ifdef xdebug 
		if ( treepartList[i].Id != i ) cout<<"ID lerde sorun var. i = " << i << endl;
#endif
/************* extra debugging ****************************************************************************/
	}

/************* time ****************************************************************************/
#ifdef showtime    		
	clock_t start = clock();
#endif    
/*****************************************************************************************/

	double potential=0;
	vektor force=0;

	Node *parentLeaf;

	//for all particles look for the local expansion and near forces contribution..
	for ( i=0; i<Nparticle;i++)
	{
		parentLeaf = treepartList[i].parent;

		//if this child is a LEAF calculate multipole expansion with partcile in this LEAF
		if ( parentLeaf->Type == LEAF )
		{
			timeTemp = clock();

			/* vec is the distance between the particle and center of root Cell*/
			vect = treepartList[i].getpos() - parentLeaf->geocenter;
			Force_C( parentLeaf->l , Mtermp , treepartList[i].getdata() , vect , &potential , &force );

			timeForceC += clock() - timeTemp;

			//we subtract the force because we calculate gravitational forces..
			treepartList[i].force -= force;
			//treepartList[i].pot += potential;

			//cout<<"ID="<<treepartList[i].Id<<" force="<<treepartList[i].force<<endl; 

			for ( k=0; k<parentLeaf->nneighbour;k++)
			{
				//calc i.particle and its near neigbour particle..
				calforcenear( i , parentLeaf->nearneighbour[k] );
				//parentLeaf is also in near neighbour of itshelf.
			}

			//if ( parentLeaf->nneighbour > 0 ) cout<<"ID="<<treepartList[i].Id<<" force="<<treepartList[i].force<<endl; 

		}//if

/************* extra debugging ****************************************************************************/
#ifdef xdebug    		 
		else 
		{
			//cout<<"Calcforce icinde cisimin parent'i LEAF degilmis.";
			cout<<i<<".cisimin parenti leaf degil=";
			treepartList[i].getpos().print();
			display( parentLeaf );			
		}
#endif
/************* extra debugging ****************************************************************************/


	}//for i.

	////////////////////////////////////////////////////////////////////////
	// now calculate the near neighbour cells particle directly.
	//////////////////////////////////////////////////////////////////////
	
/************* time ****************************************************************************/
#ifdef showtime    		
	clock_t end = clock();
	cout<<"\ncalculate force 2 take"<< ( end-start) / (float)CLOCKS_PER_SEC<<" seconds.."<< endl;
	start = clock();
#endif    
/*****************************************************************************************/

	/*when all forces calculated then calculate acceleration, velocities and positions.*/

	calcVelPos();

/************* time ****************************************************************************/
#ifdef showtime    		
	cout<<"\ntimeDirectForce = " << timeDirectForce / (float)CLOCKS_PER_SEC
		<<" -- timeForceC = " << timeForceC / (float)CLOCKS_PER_SEC <<endl;
	end = clock();
	cout<<"\ncalculate velocities and position from force take"<< ( end-start) / (float)CLOCKS_PER_SEC<<" seconds.."<< endl;
#endif    
/*****************************************************************************************/

}

/*************************************************************************
**      interaction between the source particle and the other particles  ** 
**		index is the index of array mass , Position and velocities array **
**************************************************************************/
void tree::calforcenear ( long index , Node* parentNear ) 
{
	int k;

	Node* childNode;

	if ( parentNear->Type == LEAF )
	{
		timeTemp = clock();

		//direct force calculation with all particle in childNode LEAF 
		for ( k=0; k<parentNear->Nparticle;k++)				
		{
			directforcalc( index , parentNear->plist[k]->getId() );
		}

		timeDirectForce += clock() - timeTemp;

	}//if LEAF
	//it is not necessary but control if near neighbour CELL is NODE
	else if ( parentNear->Type == NODE )
	{
////	cout<<"yakin komsusu bir NODE ve bu yakin komsunun tum cocuklari ile kuvvet hesabi cl="<<childNode->level<<" nepal="<<treepartList[index].parent->level <<endl;
		for ( k=0; k<NSUB;k++)
		{
			childNode = parentNear->child[k];
			if ( childNode!= NULL && childNode->Type != EMPTY )
				calforcenear( index , childNode );
		}
	}//else if

}

/************************************************************************************
**		calculate forces between two particle ..
**		particle particle method or Direct method...
************************************************************************************/
void tree::directforcalc( long index , long index2 )
{
	/*control that index is not equal index2*/
	if ( index == index2 ) return;
		//cout<<index<<" <--> " << index2 << endl;

	particle *part1= &treepartList[index];
	particle *part2= &treepartList[index2];

	vektor r = part1->getpos() - part2->getpos();

	double rkare = dist( part1->getpos() , part2->getpos() );

	double rkup = rkare * sqrt( rkare );

/************* extra debugging ****************************************************************************/
#ifdef xdebug    		 
	/*not Necessary but want to control collision */
	if ( rkare == 0 ) 
	{
		cout<<endl<<"fatal error distance between two particle is 0 " ;		cout<<endl;
		treepartList[index].getpos().print();		treepartList[index2].getpos().print() ;
		return;
	}
#endif
/*******************************************************************************************************/

	vektor force = part1->getforce();

	//				      M1 . M2 
	// (vektor ) F = G . --------- . (vektor) r
	//			           rkup

	force += -gravity * part1->getdata() * part2->getdata() / rkup * r;

	part1->setforce( force );

	countdirectforce++;
}


/************************************************************************************
**		calculate velocities and positions..
**		...
************************************************************************************/
void tree::calcVelPos( )
{
	/*all treeAcc is updated.*/
	for ( long i=0; i<Nparticle ; i++)
	{
		treepartList[i].calcVelPos ( dt ) ;
	}//for

	//displaytree( FMMroot );

}


/*************************************************************************
**      Find the max and min coordinate for the root cell            ** 
**      calculates the max length of one side of the root cells.      **
**  root cell coordinates vektor root , and length of one side is rsize
**************************************************************************/

void tree::findroot ( const int n , particle pos[] , double &rsize  , vektor &root )
{
   vektor maxpos , minpos , partpos;
   int i , k;

   minpos = maxpos = pos[0].getpos();            //overloaded equal operator..
   
   // find the min. and max coordinate of the system
   for ( i=1; i<n; i++ )   {
	   partpos = pos[i].getpos();
	   for ( k=0;k<NDIM;k++)
	   {
		   if ( maxpos[k] < partpos[k] ) maxpos[k] = partpos[k];
		   if ( minpos[k] > partpos[k] ) minpos[k] = partpos[k];
	   }
   }

   /* needed to calculate the length of the root cubes one side ( rsize ).  
   rsize is the max. distance between the max position and min position.
   but be careful, control all dimension one by one independent.
   */
   
   double maxdist[ NDIM ];

   for ( k=0;k< NDIM ; k++){
      maxdist[k] = maxpos[k] - minpos[k];
      //root cell coordinates.
	  root[k] = minpos[k] + maxdist[k] / 2;
   }
   
   /* find the rsize , max distance between the max. and min.*/

   rsize = maxdist[0];
   for ( k=1; k<NDIM ; k++) 
      if ( rsize < maxdist[k] ) 
         rsize = maxdist[k];

	 if ( rsize == ceil( rsize ) )
        rsize+=rsize*0.001;
     else 
        rsize = ceil( rsize );

	/* for ( k=0;k< NDIM ; k++){
      root[k] = rsize / 2;
   }*/

/************* extra debugging ****************************************************************************/
#ifdef xdebug    		 
   cout<<"\ncoordinate of center of the root cell= " << root
	   <<endl<<"and length of one side =" << rsize << endl
	   <<"region =";
   for ( k=0; k<NDIM; k++)
      cout<<"("<<root[k]-rsize/2<<","<<root[k]+rsize/2<<")";
   cout<<endl;
#endif    
/*****************************************************************************************/

   

}


/************************************************************************************
**		display root NODE or LEAF 's all information 
**
*************************************************************************************/
void tree::display( Node* root )
{
	cout.setf( ios::right );
	cout.setf( ios::fixed );
	cout.setf( ios::showpoint );
	cout.setf( ios::showpos );
	cout<< setw(4) ;
	cout<< setprecision( 2 ) ;
	
	cout<<" index="<<root->index;
	cout<<" level="<<root->level;
//	cout<<" data="<<root->data;
	cout<<" Npar="<<root->Nparticle;
	cout<<" geo="<<root->geocenter;
	cout<<" rsize="<<root->rsize;
	//cout<<"\tcms="<<root->center;
	if ( root->Type == NODE )
		cout<<"\tType=NODE"<<endl;
	else if ( root->Type == EMPTY ) 
		cout<<"\tType=EMPTY "<<endl;
	else if ( root->Type == LEAF ) 
	{
		cout<<"\tType=LEAF"<<endl;


		/*for ( unsigned int i=0; i<root->Nparticle; i++ )
		{
			
			cout<<i+1<<". particle ID= "<<root->plist[i]->getId()<<" - ";
			root->plist[i]->display();
			cout<<endl;
		}*/
	}
	//cout<<endl;

}

/************************************************************************************
**		display all tree start from root NODE toward LEAVES
**
*************************************************************************************/
void tree::displaytree( Node* root ) 
{
	int i;
	if ( root == NULL ) {
		cout<<"NULL geldi .. " <<endl;
		return;	
	}

	if ( root->Type == NODE ) 
	{
		display( root );

		for( i=0; i<NSUB; i++)				{
			if ( root->child[i] == NULL || root->child[i]->Type == EMPTY ) 
				cout<<" index = "<<i<<" level="<<root->level+1<<" is AN EMPTY or NULL NODE"<<endl;
			else //if ( root->child[i] != NULL )
			{
				displaytree( root->child[i] );
			}//if 

		}//for
	}//if
	
	// if root is a LEAF
	else if ( root->Type == LEAF ) 
		display( root );

	else if ( root->Type == EMPTY ) 
		cout<<" empty node = "<<root->index << " " << root->level;

  
	//cout<<"\nlevel of the tree = " << nlevel;

}

/************************************************************************************
**		display all tree start from root NODE toward LEAVES
**
*************************************************************************************/
void tree::displaytree2( Node* root ) 
{
	unsigned int i,k=0,j=0;
	if ( root == NULL ) {
		cout<<"NULL geldi .. " <<endl;
		return;	
	}

	for( k=0; k<NSUB;k++)
	{
		if ( root->level == j )
		{
			j++;
			display(root);
		}
	}

	for( i=0; i<NSUB; i++)				{
		if ( root->child[i] == NULL ) 
		{
				cout<<" index = "<<i<<" level="<<root->level+1<<" is AN EMPTY or NULL NODE"<<endl;
		}
		else if ( root->child[i]->Type == NODE )
		{
			//display( thisnode->child[i] );
			displaytree2( root->child[i] );
		}//if 
	}

}

